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We studied QCD with two flavors of dynamical staggered quarks at finite temperature, with a 
bare sea quark mass of about 17 MeV. We report investigations of baryon, isospin, charge and 
strangeness susceptibilities, as well as screening masses obtained from correlators of local and one- 
link separated meson operators. These were studied as functions of valence quark mass at several 
temperatures. Our results for susceptibilities deviate significantly from ideal gas values, and even 
more from the weak coupling series. We also report the first measurement of off-diagonal quark 
number susceptibilities below the transition temperature, T c , where they are the main contribution 
to charge fluctuations. We present evidence for a close connection between the susceptibilities and 
the screening masses. 
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I. INTRODUCTION 



Experiments at the Brookhavcn RHIC are now seeking to establish the formation of quark-gluon plasma. Proposals 
for experimental signatures of the plasma have to take into account its physical nature. Such information is most 
reliably obtained in lattice simulations of the hot phase of QCD. At present the phase transition temperature, T c , 
is fairly well known jl], ||, and reasonably reliable information on the equation of state has also been obtained m. 
Screening and fluctuations of conserved charges are some of the other relevant properties of the plasma. We report 
extensive new lattice results for these. 

We address two different physics questions in this paper. We are primarily interested in the question of baryon 
number and electric charge fluctuations in the plasma. These are quantities of direct experimental relevance, as has 
been realized recently [[|, ||] . We construct them from measurements of quark number susceptibilities Q . In addition, 
we study screening masses in the plasma through the study of spatial correlators, paying special attention to their 
tensor structures. This gives a new relation between the susceptibilities and screening masses which are verified by our 
lattice measurements. In particular, we explain why non-perturbative phenomena in the latter are closely connected 
with deviations from perturbation theory in the former. 

Temperatures just above T c are likely to be the most relevant region for applications in heavy ion physics and 
cosmology. Lattice measurements of the equation of state deviate significantly from the usual high temperature 
perturbative expansion. This stimulated many attempts to resum the weak coupling expansion, with varying degrees 
of success 0. Many of these techniques would have definite predictions for the susceptibilities that we measure. Since 
our measurements also deviate from the weak coupling expansion, they stand as yet another invitation to resum the 
weak coupling series. 

While the world contains six flavors of quarks, it suffices to consider two light (u and d) and a moderately heavy (s) 
quark for the physics of the QCD phase transition. Starting from a quenched theory where all quark loop effects are 
turned off, one can envisage successive approximations by which dynamical quarks are switched on in the sequence of 
their masses. Such an incremental strategy is actually forced by the cost of lattice simulations of full QCD. Earlier we 
simulated QCD in the quenched approximation, and reported our measurements of various quark susceptibilities 
In the simulations we report here, we go to the next level of approximation by using dynamical u and d quarks with 
a bare mass of about 17 MeV. Our results indicate that at temperatures above T c the change due to the inclusion of 
light quarks is small. Unquenching the heavier strange quark may thus be correspondingly less important. Even in 
the cold phase below T c , taking still lighter u and d quarks may be more crucial in future than including dynamical 
strange quarks. 

The plan of this paper is the following — technical material on the simulations, including the lattice scales and 
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details of the data taking procedure, is presented in the next section. The definitions and notation, and details of 
measurements of the screening masses and quark number susceptibilities are given in the following two sections, in 
that order. The discussion and summary in the final section is designed to facilitate use by those who are mainly 
interested in applications of our results or wish to have an overview of future directions in lattice measurements of 
susceptibilities. 

II. THE SIMULATIONS AND SCALES 
A. Run parameters 
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TABLE I: Details of Nf = 2 runs. The runs were made with a trajectory length of 1 MD time unit, a time step of 0.01 and 
a conjugate gradient stopping criterion of 10~ 5 y/V. The starred run was performed with half the MD time step but the same 
trajectory length. 
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TABLE II: Details of new quenched runs. These runs supplement earlier ones discussed in |8J. 

We have simulated QCD with two flavors of light dynamical staggered quarks (Nf = 2) with the R-algorithm || on 
N t x iVj x N z lattices with N t = 4. For investigations of quark number susceptibilities, we chose spatially symmetric 
lattices with N x — N z . On the other hand, screening correlators were easier to measure on elongated lattices with 
N z > N x . Most of our simulations consisted of generating dynamical configurations with two flavors (Nf = 2) of 
dynamical staggered quarks with the sea quark mass, m, held fixed in physical units at m/T c = 0.1 as the temperature 
was varied. A small set of runs with m/T c = 0.075 was performed to check the magnitude of sea quark mass effects. 
The run parameters and statistics collected are shown in Table [|. We emphasize that holding the quark mass fixed is 
important for physics applications, whereas the older measurements fixed m/T instead. 

After a few tuning runs which reproduced the results of previous studies |], [u|, we chose to run with parameters 
similar to those used in previous works with two flavors of staggered quarks, i.e., trajectory length of 1 molecular 
dynamics (MD) time unit, integrated in 100 steps of 0.01 MD time units each, and a conjugate gradient stopping 
criterion of \Q~ b y/V on the modulus of the remainder vector (V = N^N z N t ). It should be noted that very similar run 
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parameters are often used in Nf — 4 Hybrid Monte Carlo (HMC) runs, where a Metropolis choice ensures that the 
correct weight for a configuration is obtained even for a finite number of time steps. For Nf — 2a bias-free estimate 
of the weight is not guaranteed. In order to estimate systematic errors from this source we also made a test run with 
half the step size but for the same trajectory length. 

Part of the purpose of this work is a systematic comparison of screening masses with different numbers of dynamical 
quarks. Details of our Nr = 4 simulations are available elsewhere Jll], |l2|j . Similar details of runs in quenched QCD 
have been presented in p). However, we also generated a few sets of new quenched configurations. Details of the 
algorithm and data taking remain as in The run parameters and statistics for these quenched runs are listed in 
Table [nj. 
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FIG. 1: A diagnostic for thermalisation. The points at which the trajectory length was changed are marked by vertical arrows. 



B. Setting the scale 

Relative temperature scales, i.e., T/T c , can be deduced from previous estimates of critical couplings for various 
values of N t at the quark masses of interest to us jnj|. Thus, the temperature of a run with N t = 4 at (3 C (N[) is 
T/T c — N' t /A. Since T c is known to be 167 ± 17 MeV in physical units this allows us to deduce the value of 
lattice spacing, a, in physical units using the relation T = l/(N t a). For those values of T which cannot be obtained 
by simulations on lattices with different N t we use the two-loop beta-function of QCD to deduce the lattice spacing, 
and hence the temperature. A recent global analysis of data assures us that this is possible The values of the 
inverse lattice spacings appropriate to our measurements are given in Table |. The sea quark masses that we use are 
to = 17 ± 2 MeV (m/T c = 0.1) and 12 ± 1 MeV (m/T c = 0.075). 

Another way of setting the lattice spacing is to use previous lattice measurements of the p meson mass with the same 
quark mass and coupling |ll|. This gives us a lattice spacing 1/a = 580 ± 20 MeV at T c , and hence the dynamical sea 
quark mass comes out to be to = 14.5 ± 0.6 MeV (m/T c = 0.1) and m = 10.9 ± 0.4 (m/T c = 0.075). This consistency 
between the two estimates is, of course, not unexpected, since it is well known that T c /m p is reasonably independent 
of lattice spacing |l(J . 

The str ang e quark mass can be chosen in two different ways. One is to take the strange quark mass to be 75- 
170 MeV Q. This gives m s /T c = 0.41-1.0 [§2). Alternatively, we can use the ratio m s /m d = 17-25 @ to get 
rn s /T c = 1.7-2.5. Clearly, in the limit when the light quark masses are more realistic, the two procedures should give 
similar values. However, the heavier the strange quark mass is, the more strongly will it be subject to lattice artifacts. 
In view of this, we have chosen to work with m s /T c = 0.75, corresponding to m s — 125 MeV. A realistic ratio of the 
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strange and light quark masses will thus need future simulations with lighter u and d quark masses. Furthermore, the 
fact that the strange quark is reasonably heavy in lattice units indicates that the unqucnching effects on it would be 
small, at least near T c . 



C. Thermalisation and autocorrelations 



Due to the fact that each configuration takes a long time to generate, a major issue in any Fermion simulation is to 
identify when a run has thermalised. We monitored thermalisation through measurements of plaquettes, Wilson line, 
the chiral condensate and the number of conjugate gradient iterations needed. Since plaquette measurements are not 
very noisy, these were a good estimator of the approach to equilibrium. At temperatures above T c , the difference of 
purely spatial and mixed time-space plaquettes is nonzero, and it turns out to develop late during equilibration. This 
was one of the best criteria we observed for closeness of approach to equilibrium. However, even more efficient was 
a scatter plot of the stochastic estimator of the condensate against the number of conjugate gradient iterations 
(N cg ). As shown in Figure [l], a very definite directional movement occurs away from equilibrium, but as equilibrium is 
approached the movement becomes a random walk within a well defined area in the plane of ipip and N cg . Monitoring 
this scatter plot thus offered a possibility of tuning the algorithm to achieve fast thermalisation. 

We did this by starting each run from an ordered configuration with small trajectory lengths (t ~ 0.05 MD units). 
As the simulation proceeded, the plot of ipip against N cg began to slow. When this happened we increased the 
trajectory length in small steps until it reached unity or resulted in further motion. A post-facto check of the tuning 
was obtained after subsequent runs verified that the system continued to perform a random walk without directed 
movement. This method was developed in our earlier Nf — 4 HMC runs and worked well also for these Nf = 2 HMD 
runs. Using this technique we often managed to achieve thermalisation in 20-30 units of MD time. 

A second issue in the control of computer time is the matter of how often measurements of an observable can be 
made so that two successive measurements are effectively decorrelated. We have estimated autocorrelation times from 
a variety of gauge and fermion observables and found them to be small. Local operators such as plaquettes seem to be 
essentially decorrelated in about 2 trajectories. Even long distance observables, such as the pion screening correlator 
at the longest possible distance, are decorrelated in about 2-3 trajectories. In view of these results, we have taken 
data for screening masses on every second trajectory, and for susceptibilities every tenth trajectory. This is justified, 
a posteriori, by the fact that none of the screening masses is found to be very small. 



III. SCREENING CORRELATORS AND MASSES 



A. Symmetries and transfer matrices 



Screening involves the transfer matrix in a spatial direction of the Euclidean finite temperature lattice. Thus, 
screening correlators and masses are classified according to the symmetries of this transfer matrix, Screening masses 
are obtained from the exponential decay of the screening correlators 

c(z) = \ jjhv t Y, M W Mf ( x >y> z > t n- w 

The operators M{x,y, z,t) are made from combinations of quark and gluon fields which transform according to an 
irreducible representation (irrep) of the symmetry group, and the angular brackets denote averaging over the correct 
thermal distribution of fields. We shall also have occasion to use the meson susceptibilities 

X = ^C(z) = /^ Yl M(0)M^x,y,z,t)\. (2) 

2 \ x,y,z,t I 

Detailed discussions of the isometries of the lattice and irreps of the transfer matrix can be found for gluon operators 
in (If], 17 and staggered quark operators in |]l5| . Here we only mention the features used in this work. 



For operators made entirely out of gluon fields, "glueball" operators in the usual shorthand, the relevant symmetry 
is that of a slice of the lattice. For screening correlators such a slice includes two spatial directions and the Euclidean 
time direction. Since rotations of the spatial directions into Euclidean time are disallowed, the full cubic group, Oh, 
breaks down to the dihedral group, D\. For bilinear operators constructed from staggered quark fields, the "meson" 
operators, the symmetries of the transfer matrix are more complicated due to mixing of spin and flavor components 
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TABLE III: Representations of local and one-link separated staggered mesons that were used in this study. The table is 
extracted from jljj. Here D^x) = ^(x + fi) + <?i(x - /*)• We have chosen <q x = (-l) H+z+t , r\ v = (~l) z+t , r/ z = (-1)', r} t = 1, 
and ( x = l,£ y = (— 1)^, = ( — l) x+v , (t = (— \) x +y+ z . The spin/flavor and particle assignments at T = are standard. GiiF 
is the symmetry group of the T = staggered fermion transfer matrix and G that at T > 0. 



and staggering of the quark fields. However, this transfer matrix also carries representations of D\ as shown in [ jl5| . 
-D4 has eight one-dimensional and two two-dimensional irreps. 

In the continuum, the T = symmetry group is that of 0(3) rotations of the slice, and breaks down for T > to 
the symmetry group of the cylinder, C = 0(2) x Z2. The lattice groups are subgroups of these continuum groups. 
The real irreps of C are easily obtained. There are two one-dimensional irreps + and 0_, which come from the J z — 
components of the even and odd spin representations of 0(3). There is a countable infinity of two-dimensional real 
irreps which span the spaces of J z = ±M for any J > 0. These also carry irreps of the remaining Z2 subgroup which 
corresponds to reflections, t — * —t, in Euclidean time. 

In the continuum limit the scalar irrep At and a non-trivial irrep of D\ both collapse into the 0+ irrep of the 
cylinder group; hence these screening masses must be degenerate in this limit. The 0_ gets both the A~[ and A^ , and 
hence this pair of screening masses must also become degenerate. The four remaining one-dimensional irreps B x 2 
collapse into the J z = ±2 irrep of 0(2), and the two-dimensional irreps become the J z = ±1 of 0(2). All these 
patterns of degeneracies have been verified in the glue sector of the SU(2) pure gauge theory |PjJ. 

These seemingly esoteric group theoretic constructions actually have very simple physical applications. It has long 
been appreciated that in screening phenomena, particles of different spin may mix, and that different polarization 
components may have different dispersion relations fl§|| . These correspond, respectively, to mixing of equal J z for 
different J and the distinction between different J z for the same J. While we discuss only physics for T > and zero 
chemical potential here, it should be noted that the essential group theory lies in the inequivalence of the spatial and 
temporal directions. This is also true at finite chemical potential, and hence this group theory is also relevant to such 
future lattice computations. 



B. Results 



We have investigated the usual "local" stag gere d meson operators, i.e., those in which the quark and the anti-quark 
both land on the same lattice site (see Table III). These are familiar from their T = characters. The scalar 1 ++ , 
S, and the pseudo-scalar 1^ , PS, both belong to the At of D\. The three components of the vector 3"" ++ , V, and 
the axial- vector 3"' /H , AV, can be further reduced. These give two different copies of At, and also a Bt ■ All the 
At irreps have been measured extensively before. 

We have also made the first ever T > investigation of a class of "non-local" staggered mesons, one in which the 
quark and the anti-quark are separated by one link (see Table HI ) . These are made gauge invariant by including any 
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product of links that starts on one of these sites and ends on another. The various operators lie in the A 2 and E 
irreps of D\ . 
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TABLE IV: Screening masses in units of temperature, M/T, in various "meson" channels. A cross indicates that no data exists, 
and a dash that no stable measurement was possible. The Nf — 4 results are from runs reported in 0]. 



Chiral symmetry restoration gives the following relations between the staggered local meson correlators — 

Cps(z) = {-iyCs{z) and Cav(z) = (-l) z C v (z). (3) 

We found that these relations are satisfied for almost all the correlators to great accuracy at all T. The only small 
discrepancy is for the S/PS correlators at T = 1.5T C , where the relation is violated at the 68% confidence level (but 
not at the 95% confidence level). In order to rule out the finite MD-time step errors as the source of this oddity, we 
ran a second simulation at this coupling with half the time step. The results remained unchanged. Since our error 
estimates are stable, we performed another run with N z = 24. On this longer lattice the S and PS agreed at the 68% 
confidence level. 

For staggered Fermions, correlators are usually fitted to a form such as 

C{z) = A x cosh[Af!(z - N z /2)} + {-1) Z A 2 cosh[Af 2 (z - N z /2)]. (4) 

The correlator identities in (|J) imply that the combinations Cy ± Cav and Cy ± (— 1) z Cav project out the mass 
eigenstates in the chiral symmetric phase. We have also made single mass fits to these projections to check for 
stability of the fitted parameters. In addition we have extracted screening masses by constructing local masses with 
these projections. The local mass, m(z), is defined by the solution of 

cosh[m(z)(z - 1 - N z /2)} _ C(z - 1) 

cosh[m{z)(z ± 1 - N z /2)} ~ C(z + 1) ' ^' 

where the correlation functions are estimated by a jack-knife procedure. In order to take care of covariance between 
various measurements, we take the number of jack-knife bins to be much smaller than the number of measurements. 



In the Nf — 4 theory the B\ correlators are identically zero 11 . A statistical test of this hypothesis is made 
through the usual 

X* = - >N) (ifc " ( 6 ) 

y 

where t/j are the measurements of the correlators at different z, cr is the covariance of these measurements, and the 
hypothesis being tested is that hi = 0. With the data we have gathered in Nf = 2 QCD, we found that 

r 15.8 (1.5T C ) 

X 2 = \ 13.2 (2.0T C ) (7) 
1 24.1 (3.0T C ) 

Since N z = 16 in all these measurements, the number of degrees of freedom is 15. Thus, no alternative hypothesis 
can be supported by our data at the 99% confidence limit for T < 2T C . For T = 3T C , where there might be a 
signal, a one-mass fit by the first term of eq. (Q) gave a coefficient about 3tr away from zero, and a screening mass 
M/T = 4.5 ± 0.2. In the quenched theory, on the other hand, the situation was more like the Nf — 4 theory. At all T 
the Bt correlator was compatible with zero. The fits also gave similar agreement with a vanishing correlator. More 
detailed studies of the Bi correlator with high statistics at several temperatures above 2T C might be of interest in the 
future. 
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The correlation functions in the remaining irreps are non-trivial. The screening masses obtained are listed in Table 
IV, where we have also collected measurements made with quenched and Nf — 4 dynamical staggered Fermions on 



N t — 4 lattices. Note that the two sets of A J" correlators, the S/PS and the V/AV, give different screening masses. 
The V/AV screening mass, My, is consistent with free field theory. At all T, the screening mass is close to the 
S/PS screening mass, Ms, but due to the large errors, is also consistent with My . The E~ screening mass is harder 
to determine because the correlator is very noisy; the central value is consistent with free field theory, but has very 
large errors. 

The two parity projected local V and AV correlators are identical. This implies that not only the screening masses, 
but also the mixture of states excited by the operator are identical. In other words, the a\ and b\ screening states 
mix for T > T c . Similar equalities were found in the parity projections of the A~^ (and E~ ) correlators coming from 
the V and AV non-local meson operators, showing that the p and u> screening states also mix. The non-local S and 
PS correlators similarly bear evidence for the mixing of tt and ao- In addition, the equality of the screening mass 
determined from all four A2 correlators at our disposal also argues for a mixing of tt and the J z — component of 
the p/uj. Thus both the phenomena expected in screening correlators are seen — mixing of states of different J and 
different components of the same J having different dispersion relations. 

IV. QUARK NUMBER SUSCEPTIBILITIES 
A. Definitions 

We define the partition function 

Z = / VUe~ Sg detM(m u ,M„)det M(m d , p d ) det M(m s , p s ), (8) 



where S g is the gluon action of interest and M denotes appropriate lattice Dirac operators. The chemical potentials 
for each flavor can be combined into the singlet, triplet and octet SU(3) chemical potentials — 

Mo = M« + Pd + Ps, Pz = Pu-Pd, and p% = p u + p d - 2p s . (9) 

Our notational convention is that indices such as / and /' run over the flavors u, d and s, and indices such as i and 
j run over the SU(3) diagonal generators 0, 3 and 8. 
The quark number densities are 

where Mj. = dMf/dpf, V3 = N 3 a 3 , and T = l/N t a. Conversion to the other basis simply follows using the definitions 
in eq. (^|) and the chain rule of differentiation. Note that no = (n u + n d + n s )/3 and n 3 = (n u — n d )/2 are the baryon 
number and isospin densities, respectively. These are densities of (conserved) charges and not of quark number. The 
quark number susceptibilities are the second derivative of the free energy with respect to the chemical potentials 



dn f (T 
Xff ' ^8P^ = {V 3 



1 d 2 Z 1 dZ 1 dZ 



Z dp f dp ft Z dpj Z dpfi 



(11) 



To lighten the notation, we shall put only one subscript on the diagonal parts of \- 

We are interested in evaluating the susceptibilities for m u — m d < m s at the point pf = for all /, yielding much 
simplification. For example, each nf vanishes, a fact that we utilize as a check on our numerical evaluation. Moreover, 
the product of the single derivative terms in cq. ( |li"| ) vanishes, since each is proportional to a number density. Finally, 
since the masses are degenerate, M(m u , 0) = M(m d , 0) for each configuration of gauge links under study. 

Flavor off-diagonal susceptibilities such as 



Xud - (^r) (trM^Mi trM^M' d ) 



(12) 



are given entirely in terms of the expectation values of disconnected loops. Such quantities are numerically hard to 
compute. We discuss their evaluation later. Since M u = M d , we obtain \us = Xds with each defined by an obvious 
generalization of the formula above. Of the flavor diagonal susceptibilities we shall use 



Xs = rr 



f^j [((trAf^) 2 ) + (tr (M-'M'J - M~ 1 M' s M~ 1 ) J . (13) 
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Xu — Xd are given by an obvious generalization of this formula. Numerically, the simplest quantity to evaluate is the 
diagonal iso- vector susceptibility 



Xs = ~ (^) (tr (M-'MZ - M- l M' u M- l M' u ) ) . (14) 
Two more susceptibilities are of interest. These are the baryon number and electric charge susceptibilities, 

Xo = ^ (4%3 + Xs + 4 X«<i + 4 X«s) and x q = ^ (!0X3 + Xs + Xud - ^Xus) ■ (15) 

Note that xo is the baryon number susceptibility for three flavors of quarks. As a result, this expression differs from 
the iso-singlet quark number susceptibility for two flavors, defined in Q , both in overall normalization and by terms 
containing strangeness. 

Note that quark masses appear in two places: first in the determinant in (|^) which defines the weights for the 
averaging, and second in the trace of the Dirac operators which define the susceptibilities. The first is the dynamical 
sea quark mass m, and the second is the valence quark mass m v . In principle these can be different. We work with 
light sea quark masses for the u and d flavors, as described earlier, and use the quenched approximation for the strange 
quark, i.e., set detM s = 1. We have earlier reported measurements of these susceptibilities when all the flavors are 
quenched For staggered quarks each trace in eqs. (p2]-pT[) comes with a factor of 1/4 to compensate for flavor 
doubling. We differ from the conventions of J(| by this overall factor. 

These susceptibilities are easy to compute in free field theory i.e., for an ideal Fermi gas. As expected, nt = for 
(if = 0. All the off-diagonal susceptibilities are also zero. For X3 we obtain 

1 sin 2 p cos 2 p 

^t^x p |ro2 _|_ ^ gln 

where the spectrum of momenta is pa — (2ir/Nt)(no + 1/2), for no — 0, Nt — 1, and p v — (2tt / 'N x )n^, where 
n v = 0, • • • , N x — 1 for v ^ 0. Xs is also given by the same expression for the appropriate quark mass. Xq FT an d 
Xq FT can then be obtained using eq. (|l5|). 



B. Optimizing the measurements 



The traces required in measurements of the quark number susceptibility were estimated using a well-known stochas- 
tic method. Such stochastic estimators for the traces of N x N matrices can be constructed from ensembles of N- 
dimensional complex vectors, R, of which each component, z a , is drawn from a Gaussian distribution of unit variance. 
From the identity 



2 ' " : 



N 

n 



-M72 



2tt 



= 5 af: 



(17) 



(here d 2 z means rdrdO and the complex number z = r exp iff) we construct the estimator by multiplying both sides 
by a matrix element A a p and summing over the indices. The right hand side gives tiA. The integral on the left hand 
side has an obvious Monte Carlo estimator, giving the relation 



trA = — R\ARi - R^AR, 
2I ^ V i=i 



(18) 



where N v is the number of vectors used in making the estimate. Note that eq. (|T7]) implies that the factor 2 above is 
part of the weight in the average over the ensemble of random complex vectors. 

The stochastic estimator for (trA) 2 is a small modification of that given in M. We draw L sets of independent 
random Gaussian vectors as before and construct the estimator 



(trA) 2 = 



L(L-l) 



(RtARi) (R]AR 3 ) 



i>j=l 



(19) 



If we draw N v vectors in the evaluation of the single trace, then it is simplest to divide these into L sets of N v /L 
vectors for the estimate of each ensemble average in the sum above p3| . 
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m v /T c 


4 x 8 3 


4 x 16 3 




N v 


X3/T 2 


10 4 Xw/T 2 




X3/T 2 


W 4 X ud/T 2 


X./T 2 


0.1 


20 


0.96 (4) 


-2 (2) 


104 (6) 










40 


0.96 (3) 


-2.1 (6) 




0.98 (3) 


0.02 (2) 


106 (6) 




100 


0.96 (2) 


0.06 (4) 










0.3 


100 


0.91 (1) 


0.06 (4) 


96 (4) 


0.91 (2) 


0.02 (2) 


96 (3) 


0.5 




0.84 (1) 


0.05 (3) 


88 (3) 


0.83 (1) 


0.02 (1) 


87 (2) 


0.75 




0.748 (9) 


0.05 (3) 


78 (2) 


0.726 (8) 


0.02 (1) 


78 (2) 


1.0 




0.660 (7) 


0.04 (3) 


70 (2) 


0.639 (6) 


0.014 (8) 


70 (2) 



TABLE V: Susceptibilities for T — 1.5T C and m/T c — 0.1. On the larger lattice we used N v = 40. \t i s measured with a single 
point source by integrating the local PS correlator. 

In H the sum above was taken over all pairs, and the diagonal term was removed by subtracting a term equal to a 
stochastic estimator for tr(A 2 ). In ||] where that method was used, we showed that reasonably large values, N v w 100 
were necessary to avoid a systematic bias in Xo ~ X3- However, the estimate is facilitated by excluding the diagonal 
term from the sum above. For example, we found that the error in trA scales as 1/y/ N v , whereas that in (tvA) 2 scales 
as 1/N V , for fixed L. The dependence on L was marginal, as long N v /L > 5. 



m v /T c 


T/T c 


Xfft /T 2 


X3/T 2 


10 6 xWT 2 


X./T 2 


T/T c 


Xfft /T 2 


X3/T 2 


W 6 Xu d /T 2 


X./T 2 


0.1 


1.25 


1.1274 


0.90 (3) 


0(1) 


149 (10) 


2.00 


1.1276 


0.973 (7) 


0(1) 


84 (2) 


0.3 




1.1250 


0.80 (1) 


0.4 (8) 


119 (6) 




1.1227 


0.951 (7) 


0(1) 


82 (2) 


0.5 




1.1203 


0.704 (8) 


0.2 (7) 


99 (4) 




1.1248 


0.888 (9) 


1 (1) 


79 (2) 


0.75 




1.1112 


0.598 (6) 


0.0 (6) 


83 (2) 




1.1212 


0.849 (6) 


0.3 (8) 


75 (1) 


1.0 




1.0988 


0.511 (6) 


0.1 (5) 


72 (2) 




1.1162 


0.783 (6) 


0.4 (7) 


70 (1) 


0.1 


1.50 


1.1275 


0.96 (1) 


4(3) 


99 (4) 


3.00 


1.1277 


0.995 (4) 


0.7 (7) 


75 (1) 


0.3 




1.1259 


0.903 (7) 


3(2) 


93 (3) 




1.1273 


0.988 (4) 


0.7 (7) 


75 (1) 


0.5 




1.1226 


0.829 (6) 


3(2) 


85 (3) 




1.1264 


0.973 (4) 


0.7 (7) 


74(1) 


0.75 




1.1162 


0.733 (4) 


2(1) 


76 (2) 




1.1248 


0.946 (4) 


0.6 (6) 


72 (1) 


1.0 




1.1074 


0.646 (4) 


1.8 (9) 


69 (1) 




1.1226 


0.912 (4) 


0.6 (6) 


70 (1) 



TABLE VI: Susceptibilities measured on 4 x 12 3 lattices for Nf = 2 with m/T c = 0.1 and N v = 100. \t ls measured with a 
single point source by integrating the local PS correlator. 



m v /T c 


Xfft /T 2 


X3/T 2 


W 6 X ud/T 2 


X«/T 2 


0.075 


1.1276 


0.97 (1) 


2(2) 


106 (5) 


0.1 


1.1275 


0.97 (1) 


1 (2) 


106 (4) 


0.3 


1.1259 


0.906 (8) 


1 (1) 


98 (4) 


0.5 


1.1226 


0.828 (6) 


0(1) 


90 (3) 


0.75 


1.1162 


0.729 (5) 


0.4 (8) 


80 (2) 


1.0 


1.1074 


0.641 (4) 


0.2 (6) 


72 (2) 



TABLE VII: Susceptibilities measured at T — 1.5T C on 4 x 12 3 lattices for Nf = 2 with m/T c = 0.075 and N v = 80. 

Since the simulation time increases at least linearly in volume, we need to make a choice of volume which reduces 
computation time without introducing large artifacts into the measurements. We have made a detailed study of the 
volume dependence of susceptibility measurements at T = 1.5T C . As the data in Table M clearly show, there is no 
significant volume dependence in the values of X3 an d Xud in going from 4 x 8 3 to 4 x 16Mattices. 

In view of this, we have chosen to make the remaining measurements on 4 x 12 3 lattices. Close to T c we would 
expect that the measurements are strongly volume dependent. Based on our estimates of the screening masses we find 
that even at our smallest temperature of 1.25T C the lattice is more than 10 times larger than the longest correlation 
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m v /T c 




10 6 xWT c 2 




0.10 


0.18 (7) 


0(1) 


462 (11) 


0.30 


0.024 (9) 


0.1 (5) 


182 (3) 


0.50 


0.007 (4) 


0.1 (3) 


122 (2) 


0.75 


0.002 (2) 


0.1 (2) 


89 (1) 


1.00 


0.000 (1) 


0.0 (1) 


71.4 (8) 



TABLE VIII: Susceptibilities at T = in units of T c , measured on 12 4 lattices for N f = 2 with ra/T c =0.1 and N v = 80 



length. This assures us that finite lattice size effects are negligible at the lower end of our temperature range. On the 
other end of the scale, by the simple expedient of not going beyond 3T C , we avoid the finite volume effects such as the 
onset of spatial deconfinement. 

Finally, one must optimize the conjugate gradient stopping criterion, in which the norm of the residual vector is 
required to be smaller than e\/V. We were pleasantly surprised to find that, in thermalized test configurations, results 
for X3 and Xud changed by less than 1 part in 10 4 on changing e in the range from 0.01 to 10 -6 . Increasing e by one 
order of magnitude meant a decrease of CPU time by 20-25%. These numbers were reproduced for test configurations 
on three lattice volumes, and a variety of T. In actual measurements we chose e to be 10~ 5 for T > 1.25T C and 10~ 3 
otherwise. 



C. Results 
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FIG. 2: Xz/Xfft as a function of T/T c for m v /T c = 0.1 (boxes), 0.3 (diamonds), 0.5 (up triangles), 0.75 (down triangles) and 
1 (circles); all for m/T c = 0.1. For comparison the curves show the measurements in quenched QCD. All the data shown here 
is taken on 4 x 12 3 lattices. The shaded area denotes the mass range for strange quarks. 



Our main results on the T and m v dependence of X3 an d Xud are collected in Table VI for two flavors of dynamical 
Fermions with m/T c = 0.1. The results of a simulation at T = 1.5T C with a smaller m/T c = 0.075 are collected in 
Table VII. Comparing the two, we see no statistically significant change in X3 over this range of masses. While X3 is 
expected to be strongly dependent on m for T ~ T c in the chiral limit, our results indicate that the dependence of x 
on the sea quark mass away from the critical region is small. This suggests that the values of X3 or Xud m the chiral 
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Table VII] shows that all the susceptibilities indeed vanish 



limit for T > 1.25T C are within errors those in Table VI 
at T = as expected. 

A comparison of the dynamical and quenched theories is displayed in Figure || in terms of the ratio Xz/xfft- 
While the general trend of the data are similar in the two cases, some quantitative differences are visible. The most 
important of these is in the asymptotic value of x/xfft- For small m v the ratio is 0.85 at T/T c = 3 in the quenched 
theory, whereas it is 0.88 for Nf — 2. For larger m v , X3 is almost unchanged on unquenching the Fermions. Xs 
changes by less than 3% when the light fermions are unquenched, leading us to believe that the unquenching of the 
strange quark will not affect this quantity significantly. 

These observations are not explained in continuum perturbation theory at high temperature, which yields 



x/Xfft = 1 - 2 ( — 

V 7T 



8a 1 



6 V 7T ) 



3/2 



(20) 



when the plasmon term is resummed |19| . Here a s is the strong coupling at a scale appropriate to the temperature T. 
It is easy to see that this ratio is never less than 0.981 for Nf — (0.986 for Nf — 2). This minimum is reached when 
a s = 7r/6(6 + N f ). A recent analysis showed that TJAj^ = 1.15 ± 0.05 for Nf — and 0.49 ± 0.02 for Nf — 2 [§. 
Taking the scale for the strong coupling to be 27rT, this means that for Nf — 2 the minimum occurs at T ~ 3300A-jy^- 
(and T ~ llOA^-g- for Nf — 0). Thus in both cases, the weak coupling estimates decrease as a function of temperature 
in the range studied here, in contrast to our lattice results, which increase. It would be interesting to check whether 

a full 0(a^ 2 ) computation for x approaches the lattice results, and whether any resummations of the perturbation 
theory do better f[ p0| . 



XT O 

IT) 
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FIG. 3: Xud scales as 1/M% for T = 0.75T C in quenched QCD. The measurements of Xud and M w were both performed on the 
same configurations on 4 x 12 3 lattices ^j. 

The off-diagonal susceptibility, Xud, vanishes for T > T c and also at T — 0. As shown in Table |v|, the errors on 
Xud/T 2 for T > T c are of the order of 10 -6 , and hence this is the precision within which this quantity can be said to 
vanish for all m v . Thus, Xud and other non-diagonal parts of the flavor space susceptibilities can be totally neglected 
in constructing other susceptibilities. It is interesting to note that Xud is zero in an ideal gas, and perturbative 
contributions start at order a B . Taking the scale to be 27rT, as before, and the same values of T c /Ajj s , we might then 
expect Xud to be non-zero at the level of 0.04-0.1 (for Nf — 2) in the temperature range we studied. The substantially 
smaller, indeed vanishing, value thus accentuates the puzzle about perturbation theory [24| . 

The only place where Xud is definitely non-zero is for < T < T c . A re-analysis of our data in the quenched theory 
at T = 0.75T C l| showed that Xud differs from zero by over four standard deviations. As displayed in Figure S, Xud 
scales inversely with the square of the pion screening mass, M w , at this temperature. Since it is known that the 
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hadron masses are independent of temperature from T = to some point rather close to T c 21 , the only possible 
temperature dependence of Xud, f° r such T, would be in the value of the constant Xud,M%.. Note that the spatial 
lattice size is such that M n L > 10, and finite size effects are negligible. In order to ensure this on 4 x 12 3 lattices, we 
are constrained to use reasonably heavy quarks. The pion screening mass is then rather large. In future we plan to 
study Xud for T < T c in more detail by pushing towards lower pion mass with larger lattices, covering a bigger range 
of T and using dynamical quarks. 

D. Relation to screening masses 

The physics of a plasma lies in screening. However, it is not always clear how screening appears in different guises. 
Interestingly it can be shown that quark number susceptibility is directly related to screening, and therefore the 
non-perturbative physics in the two are likely to be the same. It is possible to present a deductive argument from 
first principles and use the data to illustrate it. Instead, we first show how the data prefers a relation between x and 
the S/PS screening mass Ms, and then give the group theoretical argument that this is not a numerical coincidence. 

The variation of X3 with m v at fixed T is shown in Figure |j. For large enough m v , there is an exponential fall of 
X3, which can be qualitatively understood as the effect of Boltzmann factors. It is less obvious that there should be a 
threshold m*(T) below which xz is almost independent of m v . The data indicates that m*(T) ~ T. Such a threshold 
cannot be derived in a weak coupling expansion. Furthermore, this threshold of constancy is exactly what prevents 
Xz/xfft from reaching unity as m v — ► 0, since Figure ^ is just another representation of the data in Figure |^. 

The only other observed non-perturbative effect in the quark sector of high-temperature QCD is in the screening 
mass in the S/PS channel, Ms- That these two effects are related is shown in the second panel of Figure ||, where we 
plot X3 against Ms/T. A simple exponential relation between them is seen. The threshold we saw before disappears 
into the relation between m v and Ms- The remaining question is why the two should be related at all. 

The answer is simple, and relies on the group theoretical classification of screening correlators and masses according 
to the symmetries of the T > transfer matrix in the spatial direction p5[ . First, from eq. (jl4|) it can be seen that 
X3 is the thermal expectation value of a product of two quark propagators sandwiching 70 U and summed over all 
distances. As a result, this quantity can be related to a susceptibility constructed from a one-link separated meson 
screening correlator. In the T = notation, this turns out to be a component of the one-link V 3 ^ For T > this 



component reduces to the irrep A 2 of the appropriate symmetry group D\ (see Table III). Other correlators which 



lie in the same irrep are the one-link separated S and PS 3" _± and the a specific component of the one-link AV 3 . 
In the continuum limit the spectrum of A 2 screening masses is degenerate with that of the , since they come from 
the same irrep of the continuum 0(2) symmetry pj]. The smallest A^ screening mass comes from the S/PS local 
propagators which are used to extract Ms- This, is the reason for the close relation between ^3 and Ms shown in 
Figure |^. 

Further numerical evidence in favor of this group theoretical argument is the similarity in the relation between 
X-it(Ms) and X3(Ms), as shown in Figure (§ The near equality of the slopes lends support to our earlier observation 
that the A 2 screening mass seems to be equal to Ms- We expect that for smaller lattice spacing or on using a Fermion 
action which restores the full flavor symmetry, the two screening masses should become equal. In this limit, the lines 
would be exactly parallel and the difference in values of x-k an d X3 at the same Ms would only reflect different operator 
overlaps with the same state. This, in fact, is a prediction which should be tested in future. 

One further piece of information fits neatly into this group theoretic framework. Xud has the same symmetry as X3 
and therefore it is also related to the S/PS correlators at all non-zero temperatures. Since it vanishes for T > T c , the 
relation is trivial there. However, for non-zero T < T c , we have observed that Xud = %/M%. This relation would be 
as mysterious as that between X3 and Ms in the absence of the argument given above. 

It is interesting to note that the preceding arguments rest entirely on the group theory of screening correlators and 
masses, i.e., on the spatial direction transfer matrix of the finite temperature system. In terms of the Euclidean time 
direction transfer matrix, the group theory is the same as at T = and there is no reason for the one-link separated 
vector temporal correlator to be related to the S/PS. Thus, the physics of quark number susceptibility is related to 
screening correlators and not to temporal correlators. 

V. DISCUSSION AND SUMMARY 
A. Event-to-event fluctuations 

Since Xq is related to charge fluctuations, it can be observed experimentally in heavy-ion collisions |1. The 
baryon susceptibility, xoi and the strange quark susceptibility, Xs, can also be used to measure fluctuations of the 
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m v /T c 




M s /T 



FIG. 4: Panel (a) shows 4x 3 /T c 2 as a function of m v /T c (for sea quark mass of 0.1T C ) at T = 1.5T C (diamonds), 2T C (circles) 
and 3T C (boxes). Panel (b) shows 4X3/T 2 (open symbols) and Xvr/IOT 2 (filled symbols) as a function of Ms/T at 2T C (circles) 
and 3T C (boxes). The lines are exponential fits. 



corresponding quantities B . Also, Xs is related by a fluctuation-dissipation theorem to the total amount of strangeness 
produced in equilibrium, and hence can be used to probe departures from equilibrium in heavy-ion collisions. Given 
the relevance of these quantities to the experimental search for the quark-gluon plasma, we have collected our results 



in an easily usable form in Table IX. 

The behavior of both Xq an d Xo above T c are essentially controlled by X3- It is clear from eq. (|l5|), that for T ^S> T c , 
where X3 — Xs Xud — Xus, we must have Xq/Xo ~ 2- Closer to T c , but still in the hot phase, if Xs is also neglected, 
we would have Xq/xo ~ 5/2. In fact, we find that this ratio decreases from 2.13 ±0.08 at 1.25T C to 2.02 ±0.01 at 3T C . 
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T/Tc 


/ FFT 
X3/X3 


/ FFT 
Xs/Xs 


/ FFT 

Xo/Xo 


/ FFT 
Xq/Xq 


1.25 
1.50 
2.00 
3.00 


0.80 (3) 
0.848 (9) 
0.863 (6) 
0.883 (4) 


0.538 (5) 
0.657 (4) 
0.757 (6) 
0.841 (4) 


0.71 (2) 
0.785 (6) 
0.828 (5) 
0.869 (3) 


0.76 (2) 
0.817 (8) 
0.846 (5) 
0.877 (3) 



TABLE IX: Susceptibilities (in units of their values for an ideal gas) obtained for degenerate dynamical u and d quarks of mass 
17 MeV and quenched strange quarks with m = 125 MeV. Susceptibilities below T c are discussed in the text. 



This shows that Xs is significant at all temperatures in the plasma, but certainly becomes closer to X3 with increasing 
T (Figure]^). In any case, charge fluctuations are twice as large as baryon number fluctuations. Interestingly, we find 
that Xq/Xs also falls from 0.95 ± 0.03 at T = 1.25T C to 0.696 ± 0.004 at T = 3T C (it is 2/3 in the T -> oo limit); 
therefore strangeness fluctuations become more significant than charge fluctuations with increasing T. 

In the low-temperature phase, < T < T c , fluctuations are dominated by the flavor off-diagonal contribution. We 
have displayed evidence (Figure ||) that Xud — l/ m w> J '- e -> fluctuations are due to the lightest particle of given flavor. 
Although the value of Xud is extremely small, the measure of fluctuation is the ratio of x an( i the entropy density [Q . 
While there are no reliable measurements of the entropy for T < T c , it is known to be small. It is presently an open 
question whether the ratio is smaller for T < T c or on the other side of the QCD phase transition. 

Assuming that the inverse relation between Xud and the pion mass generalizes to other flavor off-diagonal sus- 
ceptibilities, eq. ( |j~5|) can be used to predict that Xq/xo ~ 1/4 for T < T c . Corrections to this number are then 
given in terms of (m^/m^) 2 and (m^/m^) 2 . These might further lower the ratio by 15-20%. As a result, in the 
low-temperature phase we would obtain the hierarchy of fluctuations — xo > Xq > Xs , in total contrast to the inverted 
hierarchy Xo < Xq < Xs that we have measured above T c . 



B. Screening and susceptibilities 

In Euclidean quantum field theories at T = 0, transfer matrices in all directions are isomorphic, and one, for example, 
can focus on the time-direction transfer matrix — T = exp(— Ha) (here H is the Hamiltonian and a the lattice spacing). 
The eigenvalues of T determine the hadron spectrum. The transfer matrix and, hence, its eigenvalues do not change 
with temperature. Specifically, the symmetry of the transfer matrix remains the rotational 0(3) symmetry which is 
used to classify particle states at T = 0. 

At finite temperature or chemical potential all directions are not equivalent. As a result, T is not the same as 
any spatial direction transfer matrix. These have a totally different symmetry, C, that of the cylinder, as discussed 
extensively in the literature jlf| [l6], [l7| . If the two lowest eigenvalues belonging to the scalar of C become degenerate, 
then a phase transition occurs. In general, the eigenvalues of this transfer matrix determine screening masses in 
equilibrium. Breaking of 0(3) to C implies strange phenomena like the 'mixing' of different spins, or different states 
of the same spin taking on independent dispersion relations Jig] . While these phenomena are strange, they are not 
new — the group theory is familiar from the textbook examples of the comparison between the hydrogen atom and 
the hydrogen molecular ion . 

Lattice studies of screening correlators have hovered on the verge of this phenomenology If it has not received 
widespread attention in the past, that is merely due to the practical difficulties of measuring some of the non-trivial 



irreps of C, as we show for the first time in Table IV. The equality of parity projected correlators was used in Section 



IIIB to point out that the p and ui mix in screening, as do the ai and b\. The equality of screening masses also 



gives some evidence for the mixing of tt with the J z = componen t of th e p/u>, causing a splitting in the screening 



masses of different components of the latter. We argued in Section IVD that the observed relation between X3 and 
the screening mass M$ or between Xud and M n comes from precisely the same physics as the mixing of different spins 
at finite temperature. The same argument also implies that since fluctuations are observable, screening phenomena 
are physical. 



C. Future directions 



Several directions for future work are clear, and have been discussed in the body of the paper. Here we collect what 
seems to us the most important and fruitful possibilities. 
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One direction for future work is to examine as many of the non-local screening correlators aspossible, in order to 
gather further information on all the quantum numbers that screening correlators can come in |l5J] . As we have seen, 
this requires large lattices and immense statistics, and may well profit from the use of new noise reduction techniques 
and improved actions which do not change the symmetries (or the positivity) of the transfer matrix. 

The ratio Xz/xfft is not explained in perturbation theory. We noted this in a quenched computation and 
have verified it in the dynamical QCD computation here. The small difference (roughly 3%) between the quenched 
and dynamical results is also significant. Explanation of these results stand as invitations to those who resum 
the continuum high-temperature perturbation theory. Future lattice computations will need to push towards the 
continuum limit. Work in this direction is in progress, and will be reported soon. 

As explained before, Xud vanishes in an ideal gas of quarks, but in an interacting theory would take on a value of 
order a|. This is several orders of magnitude larger than the largest result that our measurements can tolerate. This 
observation also is an invitation to perturbation theory. 

The off-diagonal susceptibilities are non- vanishing only in the range < T < T c . This region of temperature is 
hard to study, since all the complications of T = QCD remain, and none of the simplifications of T > T c QCD start. 
On the other hand it is important input to studies of fluctuations in heavy ion collisions. We have reported a first 
quantitative observation here. Future work will push towards more realistic quark masses, taking the thermodynamic 
limit, and examining a larger range of temperatures. 
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In quenched QCD, where T c w 280 ± 10 MeV, m s /T c = 0.25-0.60. 

We thank K. Kanaya for a discussion on this estimator. 

In |^|, it is also shown that \ud starts at order a| with a small coefficient, bringing the resummed prediction closer to 
the data. 



